Viking Results, 2021-04

# Check requisite packages are installed.
packages <- c(
  "plotly", 
  "dplyr"
)
for (pkg in packages) {
  library(pkg, character.only = TRUE)
}

What do the results look like?

dirViking <- file.path(
  getwd(), "LCAB_LawMorton1996-NumericalPoolCommunityScaling"
)
dirVikingResults <- file.path(
  dirViking, "results-2021-04"
)
resultFormat <- paste0(
  "run-", 
  "%d", # Combination Number, or CombnNum.
  "-", 
  "%s", # Run Seed.
  ".RDS"
)

# Copied from LawMorton1996-NumericalPoolCommunityScaling-Calculation.R
# TODO: In the future, make this a separate file that everyone can call...
set.seed(38427042)

basal <- c(3, 10, 30, 100, 300, 1000)
consumer <- c(3, 10, 30, 100, 300, 1000) * 2
events <- (max(basal) + max(consumer)) * 2
runs <- 100

logBodySize <- c(-2, -1, -1, 1) # Morton and Law 1997 version.
parameters <- c(0.01, 10, 0.5, 0.2, 100, 0.1)

# Need to rerun seedsPrep to get the random number generation right for seedsRun
seedsPrep <- runif(2 * length(basal) * length(consumer)) * 1E8
seedsRun <- runif(runs * length(basal) * length(consumer)) * 1E8

paramFrame <- with(list(
  b = rep(basal, times = length(consumer)),
  c = rep(consumer, each = length(basal)),
  s1 = seedsPrep[1:(length(basal) * length(consumer))],
  s2 = seedsPrep[
    (length(basal) * length(consumer) + 1):(
      2 * length(basal) * length(consumer))
  ],
  sR = seedsRun
), {
  temp <- data.frame(
    CombnNum = 0,
    Basals = b,
    Consumers = c,
    SeedPool = s1,
    SeedMat = s2,
    SeedRuns = "",
    SeedRunsNum = 0,
    EndStates = I(rep(list(""), length(b))),
    EndStatesNum = 0,
    EndStateSizes = I(rep(list(""), length(b)))
  )
  for (i in 1:nrow(temp)) {
    seeds <- sR[((i - 1) * runs + 1) : (i * runs)]
    temp$SeedRuns[i] <- toString(seeds) # CSV
    temp$SeedRunsNum[i] <- length(seeds)
  }
  temp$CombnNum <- 1:nrow(temp)
  temp
})

# Note: n + 2 end states. Failure to finish, failure to obtain state, and state.
for (i in 1:nrow(paramFrame)) {
  resultsList <- list(
    "No Run" = 0,
    "No State" = 0
  )
  resultsSize <- list(
    "0" = 0
  )
  seeds <- unlist(strsplit(paramFrame$SeedRuns[i], ', '))
  for (seed in seeds) {
    fileName <- file.path(
      dirVikingResults,
      sprintf(resultFormat, paramFrame$CombnNum[i], seed)
    )
    
    if (file.exists(fileName)) {
      temp <- load(fileName)
      temp <- eval(parse(text = temp)) # Get objects.
      
      if (is.data.frame(temp)) {
        community <- toString(
          temp[[ncol(temp)]][[nrow(temp)]]
        )
        size <- toString(length(temp[[ncol(temp)]][[nrow(temp)]]))
        
        if (community == "") {
          resultsList$`No State` <- resultsList$`No State` + 1
          resultsSize$`0` <- resultsSize$`0` + 1
          
        } else if (community %in% names(resultsList)) {
          resultsList[[community]] <- resultsList[[community]] + 1
          resultsSize[[size]] <- resultsSize[[size]] + 1
          
        } else {
          resultsList[[community]] <- 1
          
          if (size %in% resultsSize) {
            resultsSize[[size]] <- resultsSize[[size]] + 1
          } else {
            resultsSize[[size]] <- 1
          }
        }
      } else {
        resultsList$`No State` <- resultsList$`No State` + 1
        resultsSize$`0` <- resultsSize$`0` + 1
      }
    } else {
      resultsList$`No Run` <- resultsList$`No Run` + 1
      resultsSize$`0` <- resultsSize$`0` + 1
    }
  }
  
  paramFrame$EndStates[[i]] <- resultsList
  paramFrame$EndStatesNum[i] <- length(resultsList) - 2 # ! No State, No Run
  paramFrame$EndStateSizes[[i]] <- resultsSize
  paramFrame$EndStateSizesNum[i] <- length(resultsSize) - 1 # ! 0
}

# X, Y, Basal and Consumer.
# Z = Sizes of the Endstates.

plotScalingData <- data.frame(
  CombnNum = rep(paramFrame$CombnNum, paramFrame$EndStatesNum),
  Basals = rep(paramFrame$Basals, paramFrame$EndStatesNum),
  Consumers = rep(paramFrame$Consumers, paramFrame$EndStatesNum)
)

# Communities
comms <- unlist(lapply(paramFrame$EndStates, names))
freqs <- unlist(paramFrame$EndStates)
freqs <- freqs[comms != "No Run" & comms != "No State"]
comms <- comms[comms != "No Run" & comms != "No State"]

plotScalingData$Communities <- comms
plotScalingData$CommunityFreq <- freqs

# Community Size
temp <- unlist(lapply(strsplit(plotScalingData$Communities, ','), length))
plotScalingData$CommunitySize <- temp

# For usage by the reader.

plotScaling <- plotly::plot_ly(
  plotScalingData,
  x = ~Basals,
  y = ~Consumers,
  z = ~CommunitySize
)

plotScaling <- plotly::add_markers(plotScaling)

plotScaling <- plotly::layout(
  plotScaling,
  scene = list(
    xaxis = list(type = "log"),
    yaxis = list(type = "log")
  )
)

plotScaling

plotScalingData

How do they compare to each other?

# > runif(1) * 1E8
# [1] 82598679
set.seed(82598679)

temp <- load(file.path(
  dirViking, 
  "LawMorton1996-NumericalPoolCommunityScaling-PoolMats.RDS"
))
mats  <- eval(parse(text = temp[1]))
pools <- eval(parse(text = temp[2]))

plotScalingData$CommunityAbund <- ""
for (r in 1:nrow(plotScalingData)) {
  plotScalingData$CommunityAbund[r] <- with(plotScalingData[r, ], toString(
    RMTRCode2::FindSteadyStateFromEstimate(
      Pool = pools[[CombnNum]], 
      InteractionMatrix = mats[[CombnNum]], 
      Community = Communities, 
      Populations = rep(100, CommunitySize), # No good guesses here
      maxRandVal = 2E4, # 2 * round of the largest value we have seen so far.
      MaxAttempts = 1E4,
      Verbose = FALSE
    )
  ))
}

Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.Failed to converge after 10000 attempts.

plotScalingData$CommunityProd <- NA
for (r in 1:nrow(plotScalingData)) {
  plotScalingData$CommunityProd[r] <- with(plotScalingData[r, ], 
    RMTRCode2::Productivity(
      Pool = pools[[CombnNum]], 
      InteractionMatrix = mats[[CombnNum]], 
      Community = Communities, 
      Populations = CommunityAbund
    )
  )
}

candidateData <- plotScalingData %>% dplyr::group_by(
  CombnNum
) %>% dplyr::mutate(
  OtherSteadyStates = dplyr::n() - 1
) %>% dplyr::filter(
  OtherSteadyStates > 0
)
candidateData %>% dplyr::select(-CommunityAbund)

islandFUN <- function(i, dat, pool, mat, dmat) {
      temp <- dat[i, ]
      RMTRCode2::IslandDynamics(
        Pool = pool,
        InteractionMatrix = mat,
        Communities = c(
          list(temp$Communities[1]),
          rep("", nrow(dmat) - 2),
          temp$Communities[2]
        ),
        Populations = c(
          list(temp$CommunityAbund[1]),
          rep("", nrow(dmat) - 2),
          list(temp$CommunityAbund[2])
        ),
        DispersalPool = 0.0001,
        DispersalIsland = dmat,
      )
    }

# For each group,
# For each pair,
# Run Island Dynamics,
# Save the result with its pairing

for (grp in unique(candidateData$CombnNum)) {
  candidateDataSubset <- candidateData %>% dplyr::filter(CombnNum == grp)
  pairingResults<- combn(
    nrow(candidateDataSubset), 2, 
    islandFUN,
    dat = candidateDataSubset, 
    pool = pools[[grp]],
    mat = mats[[grp]],
    dmat = matrix(c(0, 1, 1, 0), nrow = 2, ncol = 2),
    simplify = FALSE
  )
}

for (grp in unique(candidateData$CombnNum)) {
  candidateDataSubset <- candidateData %>% dplyr::filter(CombnNum == grp)
  pairingResults<- combn(
    nrow(candidateDataSubset), 2, 
    islandFUN,
    dat = candidateDataSubset, 
    pool = pools[[grp]],
    mat = mats[[grp]],
    dmat = matrix(c(
      0, 1, 0, # Island 2 -> 1
      1, 0, 1, # Island 1 -> 2, Island 3 -> 2
      0, 1, 0  # Island 2 -> 3
    ), nrow = 3, ncol = 3, byrow = TRUE),
    simplify = FALSE
  )
}